%function psi_h = test_stream(V,T,E,ET,TE,fe_type,dof_map,Re)
% function psi_h = test_stream(V,T,E,ET,TE,fe_type,dof_map,Re)
% stream-function vorticity form with analytic solution
%

% V T E ET TE are the mesh information
% fe_type indicate the type of finite element method, potentially maybe:

% 1 for Argyris Element 
% 2 for quintic spline the default value d = 5, r = 1 for spline space
% 3 for cubic hermit
% 4 for S_3^1(\triangle_{CT})  cubic on C-T triangulation
% 5 for S_2^1(\triangle_{PS})  cubic on P-S triangulation

Re = 40;
caseNum = 0;

% for chebyshev grid:
[V,T] = mesh_init(25,2);  % trimesh(T,V(:,1),V(:,2),'Color','blue');
[T,E,ET,TE]=mesh_edge(V,T);
bdr = sort_border(V,E,ET);  % get the boundary in counter-clockwise 

d = 5;  % the degree of spline space

% build the finite element space via D.O.F. Map
[dof_map,n_dof] = distribute_dof(T,TE,ET,d,'spline');

% and the dof constrains for C^1 spline space
[H,row_idx] = smooth_C1(dof_map,V,T,TE,ET,d);

% some auxilliary matrix for BB FEM
[cr,desc,asce,div,J,K] = auxillary_mat(d);
[Mat1,Mat2,Mat3,i_pattern,j_pattern] = ns_mat(d,desc);  % the common matrix for assembing left side

% linear part of system matrix
K = bending_mat(dof_map,V,T,d,desc); %linear part of bending matrix

% setup the parameters for the equation solver
epsilon = 5e-4; 
Max_linear_iter = 100; linearTol = 1e-5; 
Max_newton_iter = 12; newtonTol = 1e-8;

% 1. Calculate right hand sides now
b = rhs_bb(V,T,d,dof_map,desc,Re,caseNum);

% 2. Boundary condition, this time, it is dominanted
%[B,G] = get_navstk_bc(V,T,TE,ET,d,dof_map,cr,bdr,Re,caseNum);
[B,G1] = inflow_bc(V,T,TE,ET,d,dof_map,cr,bdr,desc,Re,caseNum);

% solve the Stokes problem for initial value;
p = size(H,1); ZERO = zeros(p,1); L = [H;B]; G = [ZERO;G1];

% 
c0 = lagrange22(K/Re,b,L,G,epsilon,Max_linear_iter,linearTol);

% solve the NavStk with the Newton method
err = 50; iter = 1; dim = n_dof;
fprintf('\nTrying Renold = %f, system sized %d x %d.\n',Re,dim,dim);
while err > newtonTol && iter <= Max_newton_iter
 
    [Q,QP] = NS_Newton(V,T,d,dof_map,c0,Mat1,Mat2,Mat3,i_pattern,j_pattern); 
    DF = K/Re + QP;
    F = Q + b;
    c = lagrange22(DF,F,L,G,epsilon,Max_linear_iter,linearTol);

    diff = c-c0; err = norm(diff,inf);
    c0 = c;  % prepare old value for next newton iteration
    fprintf('Iteration %d error : %e \n',iter,err);
    iter = iter + 1;

end

cvg = 0;
if err <= newtonTol   % convergent , then plot the solution size(X)),L);
    cvg = 1;
end

xmin = min(V(:,1)); xmax = max(V(:,1));
ymin = min(V(:,2)); ymax = max(V(:,2));
x = linspace(xmin,xmax,100); y= linspace(ymin,ymax,100);

[vtx,tri,ph] = heval(dof_map,V,T,c,d);
[p,px,py] = psi_h(vtx(:,1),vtx(:,2),Re,0);
trisurf(tri,vtx(:,1),vtx(:,2),ph-p);


%[X,Y] = meshgrid(x,y);
%Z = geval(dof_map,V,T,d,c,X(:),Y(:));
%Zh = psi(X,Y,Re);
%surf(X,Y,reshape(Z,size(X)));%
%norm(reshape(Z,size(X))-Zh)
% shading interp;  % min(min(Z))-max(max(Z))

% L1 = [-0.11,-0.9,-0.7,-0.5,-0.35,-0.1,-1e-5,1e-5,1e-3,1e-2];
% % L=[-1e-10 -1e-7 -1e-5 -1e-4 -0.01 -0.03 -0.05 -0.07 -0.09 -0.1 -0.11 -0.115 ...
% %                 -0.1175 1e-8 1e-7 1e-6 1e-5 5e-5 1e-4 2.5e-4 5e-4 1e-3 1.5e-3 3e-3 0 0.5  ...
% %                 1 2 3 4 5 -0.5 -1 -2 -3];

% figure;contour(X,Y,reshape(Z,size(X)),40);
% hold on;contour(X,Y,reshape(Z,size(X)),L1);
% pos = find(Z==min(min(Z)));
% plot(X(pos),Y(pos),'r*'); hold off;